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ABSTRACT 


The  approximation  of  definite  integrals  using  Monte  Carlo  simulations  is  the 
focus  of  the  work  presented  here.  The  general  methodology  of  estimation  by 
sampling  is  introduced,  and  is  applied  to  the  approximation  of  two  special 
functions  of  mathematics:  the  Gamma  and  Beta  functions.  A  significant  ap¬ 
plication,  in  the  context  of  radar  detection  theory,  is  based  upon  the  work  of 
[Shnidman  1998].  The  latter  considers  problems  associated  with  the  optimal 
choice  of  binary  integration  parameters.  We  apply  the  techniques  of  Monte 
Carlo  simulation  to  estimate  binary  integration  detection  probabilities. 
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Approximation  of  Integrals  via  Monte  Carlo  Methods,  with 
an  Application  to  Calculating  Radar  Detection  Probabilities 


EXECUTIVE  SUMMARY 


The  performance  analysis  of  a  radar  detection  scheme  requires  estimation  of  probabilities 
of  false  alarm  and  detection,  under  various  clutter  scenarios.  These  probabilities,  which 
often  appear  as  definite  integrals,  are  frequently  analytically  difficult  to  evaluate.  Hence, 
numerical  approximation  schemes  are  employed.  Monte  Carlo  estimators  use  statistical 
simulation  to  evaluate  such  integrals.  As  with  any  approximation  scheme,  there  are  limi¬ 
tations  and  drawbacks  in  its  application.  One  of  the  major  difficulties  with  Monte  Carlo 
estimators  is  that  very  large  sample  sizes  may  be  required,  in  order  to  achieve  a  reasonable 
estimate.  This  is  especially  true  in  the  context  of  estimating  probabilities  of  rare  events, 
such  as  radar  false  alarms. 

The  purpose  of  this  report  is  to  examine  the  Monte  Carlo  estimation  of  integrals  in  general. 
After  formulating  the  scheme,  applications  to  the  evaluation  of  two  special  functions  are 
considered.  The  success  of  an  estimator  will  be  decided  on  its  performance  in  terms  of 
providing  a  reasonable  estimate  for  the  smallest  sample  size  possible. 

The  major  application  in  this  report  will  be  to  obtain  estimates  for  a  detection  probability 
in  a  binary  integration  context.  Under  the  assumption  of  a  Swerling  target  model,  an 
expression  for  the  binary  integrated  probability  of  detection  is  obtained  in  Shnidman’s 
1998  paper  entitled  Binary  integration  for  Swerling  target  fluctuations  (IEEE  Transactions 
on  Aerospace  and  Electronic  Systems,  Volume  34,  pp.  1043-1053).  We  apply  Monte  Carlo 
simulations,  together  with  some  functional  approximations,  to  estimate  this  probability 
for  Swerling  1  and  3  target  models. 
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1  Introduction 


It  is  a  common  occurrence,  in  the  study  of  radar  detection  performance,  to  find  it  ana¬ 
lytically  intractible  to  construct  a  closed  form  expression  for  probabilities  of  interest.  The 
two  key  performance  measures,  the  probability  of  false  alarm  and  probability  of  detection, 
used  in  the  analysis  of  Constant  False  Alarm  Rate  (CFAR)  detectors,  typically  involve 
integrals  that  cannot  be  readily  evaluated.  Monte  Carlo  methods  are  thus  often  used, 
and  in  the  context  of  false  alarm  probabilities,  much  work  has  been  generated  on  the 
construction  of  suitable  suboptimal  biasing  densities  (see  [Weinberg  2004]  and  references 
contained  therein).  In  a  CFAR  system,  the  false  alarm  probability  is  set  to  a  very  small 
number,  and  consequently  Monte  Carlo  estimators  of  this  may  need  a  very  large  sample 
size  to  achieve  a  nonzero  estimate. 

In  the  context  of  CFAR  detection,  we  are  testing  whether  a  test  observation  xq  represents 
a  target  or  not.  This  decision  is  based  on  whether  this  observation  exceeds  a  weighted 
“average”  measure  of  the  clutter  level.  Mathematically,  if  So  is  the  random  variable 
representing  the  test  observation,  and  Si,  S2, . . . ,  S„  are  clutter  statistics,  we  declare  a 
detection  if  So  >  ^^(Si,  S2, . . . ,  S „),  where  r  is  the  threshold,  and  is  a  function  which 
measures  the  clutter  level.  The  probability  of  false  alarm  and  detection  can  be  written  in 
the  form 


P 


rJr 


■  J  G(xq,xi,...,  xn)£(x0,xi, . . . ,  xn)dxodxi  •■■dxri 


Here  C(^o?  #1 ,  •  •  •  5  xn)  —  I[xo  >  ~'ip(xi,x2i  •  •  • ,  xn)]i  where  I  is  defined  by 


l[x  e  A] 


1  if  x  £  A; 

0  otherwise. 


(i) 


The  term  involving  £  is  the  joint  density  of  the  cell  under  test  and  the  clutter  statistics. 
Whether  the  integral  (1)  is  for  a  detection  or  false  alarm  probability  will  depend  on  the 
distribution  of  the  cell  under  test  statistic.  In  either  direction,  what  becomes  apparent  is 
that  this  integral  will  often  be  quite  difficult  to  evaluate  analytically.  Hence  a  numerical 
approximation  scheme  is  required.  Due  to  the  presence  of  a  density  in  the  integral  (1), 
Monte  Carlo  estimation  seems  to  be  a  natural  choice. 

The  objective  of  this  report  is  to  illustrate  the  application  of  Monte  Carlo  methods  to  the 
more  general  problem  of  integral  evaluation.  The  ideas  of  changing  a  simulation  distribu¬ 
tion,  also  known  as  Importance  Sampling ,  will  be  illustrated  in  this  context.  Monte  Carlo 
techniques  will  be  illustrated  in  the  application  to  evaluation  of  some  special  functions  that 
arise  in  mathematics.  A  specific  radar  related  application  appears  in  the  evaluation  of  a 
detection  probability  integral.  The  latter  arises  in  the  context  of  binary  integration  with 
Swerling  target  models,  and  appears  in  [Shnidman  1998].  A  Monte  Carlo  scheme  is  used, 
as  well  as  some  other  approximations,  to  construct  estimates  for  the  binary  integrated 
detection  probability  for  two  Swerling  models. 

We  begin  by  introducing  the  basic  ideas  of  Monte  Carlo  methods. 
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2  Monte  Carlo  Techniques 


The  application  of  simulation  methods  to  the  estimation  of  difficult  integrals  began  with 
the  work  of  [Kahn  1950]  and  others  working  in  nuclear  physics  in  the  1940s.  An  early 
radar  application  to  estimation  of  false  alarm  probabilities  is  [Mitchell  1981].  The  basis 
for  Monte  Carlo  techniques  is  the  Strong  Law  of  Large  Numbers  (SLLN)  (see  [Billingsley 
1986]).  This  states  that  a  series  of  independent  and  identically  distributed  (HD)  random 
variables,  normalised  by  the  number  of  terms  in  the  series,  will  converge  to  the  mean  of 
any  one  of  the  terms  in  the  series.  Mathematically,  this  implies  that  if  Hi,  E2, . . . ,  is  a 
sequence  of  IID  random  variables  with  finite  mean  i?[E],  then 


m 

^3 

lim  ^ —  =  i?[E] , 

m — >00  yyi 

except  on  a  set  of  probability  measure  zero.  This  suggests  that,  for  sufficiently  large  m, 
the  average  of  the  random  variables  in  (2)  can  be  approximated  by  its  mean.  Where  this 
applies,  in  the  context  of  interest,  is  that  it  enables  the  evaluation  of  integrals  through 
simulation. 


Consider  the  integral  I  =  f^u(x)^(x)dx^  where  £  is  a  density  on  $7,  and  a;  is  a  function. 
This  integral  is  a  statistical  mean  of  a  random  variable  E  on  £1  with  density  £.  Hence  we 
can  write  I  —  £/[o;(E)].  Now  if  Ei,E2,...,Em  is  an  IID  sequence  of  random  variables, 
then  the  SLLN  (2)  implies  that 


£  w(sj) 

lim  - — - - =  l?[w;(E)l  =  I,  (3) 

m— >oo  777, 

except  on  a  set  of  probability  measure  zero.  Hence,  by  applying  (3),  we  can  deduce  that 

m 

I  =  [  u(x)£(x)dx  &  - ,  (4) 

Ja  rn 

where  the  sequence  22, . . . ,  consists  of  realisations  of  the  variables  Si,  E2, . . . ,  Em. 
The  result  in  (4)  implies  that  the  integral  I  can  be  estimated  by  generating  a  series 
of  realisations  of  a  random  variable  with  density  £,  and  evaluating  the  average  of  the 
function  u  over  these  realisations.  This  is  a  computationally  simple  exercise  in  theory. 
As  remarked  previously,  an  underlying  problem  with  Monte  Carlo  methods  is  that  it  may 
require  a  very  large  sample  size  to  achieve  a  reasonable  estimate.  Changing  the  biasing 
density  can  sometimes  rectify  this,  and  this  will  be  considered  in  the  discussion  to  follow. 

In  cases  where  we  have  a  definite  integral  of  a  function  that  is  not  a  density  on  the  integral’s 
domain,  it  is  still  possible  to  apply  Monte  Carlo  methods.  To  illustrate  this,  suppose  I  is 
a  definite  integral  of  the  form 

I  =  [  C (x)dx,  (5) 

**  sA. 
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where  £  :  A  — ►  R  is  a  continuous  nonnegative  real  valued  function  on  the  interval  A  = 
[a,  b\.  We  do  not  assume  that  the  boundaries  of  this  interval  are  finite,  so  that  we  allow 
for  integrals  on  infinite  domains.  We  do  assume,  however,  that  the  integral  exists,  in  a 
Riemann  or  Lebesgue  sense.  We  would  like  to  apply  the  SLLN  to  (5),  in  order  to  apply  a 
Monte  Carlo  approximation.  We  can  construct  a  probability  density  function  £  on  A,  and 
modify  the  integral  to 

I  =  /  jr\i{x)dx  =  JEw(H),  (6) 

Ja  €(.x) 

where  u(x)  =  and  the  expectation  in  (6)  is  with  respect  to  a  random  variable  H  with 
density  £.  We  require  this  density  to  be  nonzero  so  that  u(x)  is  well  defined.  The  existence 
of  such  densities  can  be  shown  by  considering  the  four  types  of  integral  domains.  If  the 
integral’s  domain  is  an  interval  of  the  form  [a,  6],  where  both  a  and  b  are  finite,  then  one 
can  choose  a  uniform  distribution  over  this  domain.  In  the  case  of  an  integral  with  domain 
[a,  oo)  or  (— oo,  ft],  where  both  a  and  b  are  finite,  an  Exponential  distribution  can  be  used. 
The  final  case,  where  the  integral  is  over  the  whole  real  line,  a  Gaussian  distribution  can 
be  used.  This  procedure  is  often  referred  to  as  Importance  Sampling.  [Weinberg  2004] 
contains  a  detailed  list  of  Importance  Sampling  references. 


An  application  of  the  SLLN  to  (6)  results  in  the  Monte  Carlo  estimator 


In 


N 


N  ^  K 


3= 1 


(7) 


where  each  H j  =  S.  In  general,  the  expression  $  =  means  that  the  distributions 
of  random  variables  <&  and  are  equal,  so  that  for  every  set  A  in  a  common  domain, 
P(4>  E  A)  =  P(\?  E  *4).  Estimator  (7)  is  an  unbiased  estimator  of  /,  since 


IE Jjv  =  IEcj(S)  =  [  u(x)£(x)dx  -  I.  (8) 

Ja 

Thus  estimates  are  centred  on  the  integral  being  approximated.  The  variance  of  (7)  can 
be  shown  to  be 


WIN 


I  (lEMS)2)  -  (IE(u,(E))2) 


(9) 


We  would  like  to  have  an  estimator  whose  variance  is  as  small  as  possible.  Notice  that 
with  the  choice  of  £  (x)  =  on  A,  the  variance  (9)  reduces  to  zero.  This  biasing  density 
is  known  as  the  optimal  solution ,  but  is  of  no  practical  use  because  it  depends  on  the 
quantity  being  estimated.  Many  authors  have  used  knowledge  of  the  optimal  solution  to 
construct  a  suboptimal  biasing  density,  with  varying  degrees  of  success  (see  [Gerlack  1999] 
and  [Orsak  and  Aazhang  1989]).  Its  form  suggests  that  a  suitable  biasing  density  should 
be  concentrated  on  the  integral’s  domain,  and  in  some  sense  proportional  to  the  integrand. 
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Thus,  if  we  are  presented  with  an  integral  of  the  form  (5),  there  are  three  ways  to  proceed 
in  terms  of  the  Monte  Carlo  approach.  If  integrand  £{x)  contains  a  density  on  the  integral’s 
domain,  we  can  use  this  to  simulate.  If  this  is  not  the  case,  then  we  can  insert  a  biasing 
density,  and  alter  the  integrand.  The  other  possibility  is  that  we  can  change  biasing 
densities  even  if  (  does  contain  a  density. 


We  now  turn  to  the  issue  of  a  suitable  choice  of  biasing  density.  There  are  many  different 
distributions  from  which  one  can  select  a  biasing  distribution.  The  Weibull  family  of 
distributions,  yV(a,/?),  with  nonnegative  parameters  a  and  /?,  has  probability  density 
function 


Oi  ( X 

p  \p 


OL—1 


? 


and  can  be  simulated  via 


=/?(-log(ii))“  =  W(«,0), 

where  <J>y\;(:r)  is  the  cumulative  distribution  function,  and  R  =  R( 0, 1)  is  a  continuous 
uniform  distribution  on  the  unit  interval.  There  are  a  number  of  important  special  cases 
of  the  Weibull  distribution.  Observe  that  W(1 ,/?)  is  an  Exponential  distribution,  while 
W(2,  (3)  is  Rayleigh. 


The  Gamma  family  of  distributions,  <?(r,  /?),  has  density 


_ x_ 

e  v 


Pi  (r) 


where  7(7*)  is  a  normalising  constant,  called  the  Gamma  function,  and  we  also  assume 
parameters  r  and  j3  are  nonnegative.  It  is  analytically  impossible  to  write  down  the 
inverse  of  the  cumulative  distribution  function  of  the  Gamma  distribution.  However, 
it  is  possible  to  simulate  from  this  distribution,  using  the  fact  that  if  H  is  a  random 
variable  with  the  Gamma  distribution  with  parameters  r  and  /?,  then  S  =  H1  +  S2  + 
•  •  •  +  Hr,  where  each  5 j  is  an  IID  Exponential  random  variable  with  parameter  f3.  Thus 
realisations  of  a  Gamma  distribution  can  be  obtained  by  summing  independent  realisations 
from  Exponential  distributions. 


The  success  of  a  biasing  density  is  largely  application  dependent.  In  the  next  section 
we  will  investigate  the  result  of  applying  different  biasing  densities  to  the  Monte  Carlo 
estimation  of  some  special  functions. 
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3  Monte  Carlo  Approximations  of  some  Special 

Functions 


In  order  to  illustrate  the  application  of  Monte  Carlo  simulation  to  the  evaluation  of  inte¬ 
grals,  we  examine  two  special  functions.  The  Gamma  and  Beta  functions  arise  in  a  number 
of  contexts  in  mathematics  and  statistics.  In  applied  mathematics,  these  functions  appear 
in  the  theory  of  Bessel  functions  [Bowman  1958].  They  arise  as  normalising  constants  in 
the  definition  of  certain  distributions  in  the  theory  of  probability  [Ross  1983].  The  Gamma 
function  7(71)  :  R+  — >•  R  is  defined  by  the  integral 

poo 

7  (n)  =  /  xn~] e~x  dx,  (10) 

Jo 

and  the  Beta  function  /?(ra,  n)  :  R+  x  R+  — *•  R  is 

f3(m,n)  =  i  —x)n~1dx.  (11) 

Jo 

For  nonnegative  integral  n,  it  can  be  shown  that  7 (n  +  1)  =  n!.  The  Beta  function  can  be 
decomposed  into  an  expression  involving  Gamma  functions.  Specifically,  it  can  be  shown 
that  /3(m,  n)  =  .  Except  for  some  simple  cases,  it  is  necessary  to  use  numerical 

methods  to  evaluate  these  integrals  in  practice.  We  will  apply  Monte  Carlo  techniques  to 
both  (10)  and  (11).  In  the  case  of  (10),  the  exponential  term  £(#)  :=  e~x  is  a  density  on 
the  interval  [0,  00).  Hence,  an  appropriate  Monte  Carlo  estimator  of  (10)  is 

1  N 

7i(n,lV)  -  —  »  (12) 

3= 1 

where  each  5  is  a  random  variable  generated  from  a  distribution  with  density  £.  In  this 
case,  since  £  is  the  density  of  an  Exponential  random  variable  with  parameter  1,  we  can 
use  5  =  —  log(i?)  to  simulate  the  distribution,  where  R  is  a  uniform  density  on  the  interval 
[0,1]. 


We  can  construct  an  alternative  estimator  by  simulating  from  a  distribution  with  density 
from  the  Weibull  class.  If  £  is  the  Weibull  density  £(x)  =  e  ^ '  ,  then  an 

estimator  based  on  this  is 


72  (n,  N) 


1  N 
-T 

N  ^ 


3= 1 


R°i  N  _C7  . 

P  V-  ~n-OL  ^ 
aN  ^  ^ 

3= 1 


where  each 


~j  =  j3(—  log  Rj ) « ,  with  Rj  =  R[ 0, 1]. 


(13) 


Figures  1-3  in  Appendix  A  contain  simulations  to  estimate  the  Gamma  function,  using 
both  the  estimators  (12)  and  (13).  The  Matlab  code  used  to  generate  the  estimates  can  be 
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found  in  Appendix  B.  Each  plot  compares  these  two  estimators  to  the  exact  value  for  the 
Gamma  function,  or  a  numerical  approximation  based  upon  an  inbuilt  Matlab  function. 
Figures  1  and  2  are  based  upon  a  sample  size  of  N  =  10,  000,  while  that  of  Figure  3  is 
based  upon  a  sample  of  size  of  N  =  100,  000.  The  Weibull  estimator  (13)  was  chosen  with 
parameters  a  =  1  and  (3  =  10.  It  was  found  empirically  that  other  values  for  a  did  not 
provide  good  estimates.  The  simulations  show  there  is  reasonably  good  convergence  to 
the  exact  value  for  approximately  100,000  simulation  runs. 

The  Beta  function  can  also  be  estimated  using  Monte  Carlo  estimators.  The  difference  is 
that  the  domain  of  the  integral  is  the  unit  interval  [0,1].  A  suitable  density  on  this  interval 
is  the  standard  uniform  one,  which  is  £(x)  =  1,  for  x  E  [0, 1].  In  this  case,  the  Monte 
Carlo  estimator  is 

_  i  N 

n,  N)  =  —  Y,  (!  “  >  (14) 

3=1 

where  in  this  case,  each  S j  is  a  random  variable  with  the  standard  uniform  density  on 
[0,1].  There  are  also  some  other  natural  choices  for  biasing  densities  for  the  Beta  function. 
We  can  choose  £(#)  =  ra:rm_1,  and  introduce  a  scaling  factor  of  m  to  the  Beta  function’s 
integral.  Hence  the  estimator  is 

_  1  * 

/32(m,n,N)  = \  (15) 

3= 1 

where  5  is  a  random  variable  with  the  density  £.  This  can  be  simulated  from  the  fact 
that  R™  =  S,  where  R  =  R[ 0,1]  is  a  uniform  distribution  on  the  unit  interval.  Note  that 
a  biasing  density  can  also  be  based  on  the  second  term  in  the  integrand.  Specifically,  we 
could  make  the  choice  of  £(#)  =  n(l  —  x)n~ which  is  also  a  density  on  the  unit  interval. 

As  remarked  in  Section  2,  it  is  possible  to  simulate  from  any  distribution  we  can  define 
on  the  integral’s  domain.  In  the  current  context,  we  could  insert  a  modified  Exponential 
distribution  into  the  Beta  integral,  and  use  this  for  simulation.  To  illustrate,  note  that 
the  function  £(#)  =  f  is  a  density  on  [0, 1],  and  a  random  variable  S  with  this  as  its 

density  can  be  simulated  by  —  log(l  —  (1  —  e~1)R)  =  S,  where  as  before  R  is  the  standard 
uniform  distribution  on  the  unit  interval.  This  distribution  is  referred  to  as  a  truncated 
Exponential  distribution.  Thus  a  third  estimator  of  the  Beta  function  is 

fo(m,n,N)  =  — (16) 

3=1 

where  Ej  is  generated  from  the  truncated  Exponential  distribution. 

Figures  4-6  in  Appendix  A  contain  a  number  of  simulations  of  the  three  estimators  (14)- 
(16).  For  given  values  of  m  and  n,  these  estimators  are  compared  to  the  exact  result. 
In  contrast  to  the  estimates  for  the  Gamma  function,  there  is  rapid  convergence  in  this 
case.  After  only  100  simulations,  the  estimators  are  giving  very  good  estimates  of  the  Beta 
function. 
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4  Monte  Carlo  Approximations  of  Detection 

Probabilities 

We  now  turn  to  the  problem  of  estimating  radar  detection  probabilities  using  Monte  Carlo 
methods.  The  detection  probability  under  consideration  arises  in  the  context  of  binary 
integration  for  Swerling  target  fluctuations  [Shnidman  1998].  As  an  alternative  to  coherent 
integration  of  N  pulses  [Levanon  1988],  binary  integration  determines  how  many  single 
pulse  exceedences  of  the  threshold  occur.  A  target  detection  is  declared  if  there  are  at 
least  M  such  exceedences,  for  a  prescribed  value  of  M,  within  N  observations.  A  key 
issue  is  to  determine  the  optimal  value  of  M,  which  is  the  focus  of  [Shnidman  1998].  We 
will  instead  assume  such  an  optimal  value  has  been  determined,  and  will  investigate  the 
corresponding  detection  probability  integrals. 

The  following  is  taken  from  [Shnidman  1998],  with  slight  modification  of  notation.  Define 
the  cumulative  distribution  function  of  Binomial  probabilities  between  two  integers  N  and 
M,  with  1  <  M  <  Af,  to  be 

E(N,  M,p)  =  Y,  (fy  (i  -  p)N~kpk,  (17) 

for  some  0  <  p  <  1.  The  binary  false  alarm  probability  is  Pfa  —  -E^AT,  M,pi),  where 
pi  is  the  single  pulse  false  alarm  probability.  The  single  pulse  normalised  threshold  is 
r  =  —  log(pi).  We  define  5  to  be  the  target  normalised  signal  to  noise  ratio  (SNR),  which 
under  Swerling  models,  we  assume  has  a  Gamma  distribution  with  parameters  r  =  l  and 
(3  =  y.  The  parameter  l  is  a  fluctuation  parameter,  and  £o  is  the  average  normalised 
SNR.  The  density  of  S  is  thus 

(18) 

We  let  pc(t ,  r)  be  the  single  pulse  probability  of  detection,  for  a  constant  target  with  single 
pulse  normalised  SNR  level  t.  Then 


where  Iq  is  the  modified  Bessel  function  of  order  zero  [Bowman  1958].  For  the  four 
Swerling  target  models  considered  in  [Shnidman  1998],  the  binary  integrated  probability 
of  detection  turns  out  to  be 

r  oo 

PD=  /  (20) 

Jo 

where  the  form  of  p(t,r)  depends  on  the  Swerling  case.  As  explained  in  [Shnidman  1998], 
for  Swerling  1  and  3,  p(t,  r)  =pc(t,r).  In  the  case  of  Swerling  2, 

roo  _  t 

P(t,r)=  pc(p,T)Z(p\Zo,l  =  l)dp  =  e  x+«o ,  (21) 

Jo 
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where  the  latter  equality  can  be  demonstrated  analytically.  It  is  important  to  note  that 
the  fluctuation  parameter  l  is  taken  to  be  1,  for  Swerling  1  and  2  target  models,  and  is 
assumed  to  be  2  for  Swerling  cases  3  and  4.  In  the  case  of  Swerling  4, 

roo 

p(t,r)=  pc(p,T)£(p\£o,l  =  2)dp,  (22) 

J  0 

the  difference  between  this  and  (21)  being  the  difference  in  fluctuation  parameters.  When 
applied  to  (20),  the  Swerling  2  and  4  expressions  for  p(t,  r)  do  not  depend  on  t,  and  so  (20) 
reduces  to  the  cumulative  sum  of  binomial  probabilities  (17).  Hence  we  do  not  examine 
these  cases  further,  since  Monte  Carlo  simulations  are  not  required.  We  would  like  to 
obtain  estimates  of  (20),  using  a  Monte  Carlo  approximations,  for  Swerling  1  and  3  target 
models. 

The  presence  of  the  Gamma  density  in  (20)  suggests  that  this  could  be  used  as  a  simulation 
density.  Hence,  a  Monte  Carlo  estimator  for  the  detection  probability  (20)  is 

1  K 

PD(N,  M,K,t)  =  -J2  E(N’  r)),  (23) 

3= 1 

where  each  H j  is  a  Gamma  random  variable  with  density  given  by  (18).  The  Gamma 
variables  can  be  simulated  by  adding  r  =  l  independent  realisations  of  Exponential  random 
variables  with  parameter  (3  =  y. 

We  now  need  to  approximate  pc(t,r).  There  are  two  possiblities.  Firstly,  from  [Shnidman 
1995]  and  [Weinberg  and  Kyprianou  2005],  we  have 


Pc(t,r) 


1  +  i(l 


A 

3! 


(l  +  r  + 


2! 


+  .  .  .  , 


(24) 


which  yields  the  approximation 

(  r2\2 12 

log(pc(t,r))  fa  -T  +  rt-  (  r  +  —  j  — .  (25) 

Note  that  this  quadratic  expression  is  always  negative  which  is  important  since  pc(t,r)  is 
a  probability,  and  any  estimate  of  it  must  also  be  a  probability. 

We  would  expect  the  approximation  (25)  to  not  work  well  except  when  t  and  r  are  both 
small. 
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Secondly,  from  [Weinberg  and  Kyprianou  2005],  we  have 


Pc(t,r)  =  P(X2(t)  <Xi(t)), 


(26) 


where  (^i(t),  ^(r))  =  (Po(£),  P°(r))  are  independent  Poisson  random  variables.  (Note 
that  Po(A)  indicates  a  Poisson  distribution  with  mean  value  A). 

A  Monte  Carlo  estimator  of  (19)  can  be  based  on  (26),  by  simulating  from  the  distribution 
of  Hi  and  evaluating  cumulative  probabilities  of  ^2*  Specifically,  one  can  use  the  estimator 

1  H 

Pc(t,  t)  =  —J2  P(*2(t)  <  £,),  (27) 

3= 1 

where  each  is  generated  from  a  random  variable  with  a  Po(t)  distribution.  This  can  be 
used  in  conjunction  with  the  estimator  (23)  to  estimate  the  detection  probability  (20).  It 
was  found  that  an  estimator  based  on  (27)  converged  faster,  and  for  less  simulation  runs, 
than  an  estimator  based  directly  on  the  integral  (19).  Using  (27)  also  provided  better  sim¬ 
ulation  estimates  for  the  probability  (20),  rather  than  using  the  quadratic  approximation 
(25). 

Figure  7  in  Appendix  A  is  a  simulation  of  (23),  in  the  case  where  TV  =  5,  M  —  3, 
l  —  2,  £o  —  1  and  T  =  0-4.  The  estimator’s  Matlab  code  can  be  found  in  Appendix  B. 
The  quadratic  approximation  (25)  was  used  to  estimate  the  probability  (19).  The  ^th 
simulation  uses  K  =  10J  iterations  in  the  estimator  (23).  As  can  be  observed,  the  Monte 
Carlo  estimator  (23)  begins  to  settle  down  from  the  third  simulation,  which  corresponds 
toif=  1000. 

Figures  8-11  in  Appendix  A  contain  a  number  of  simulations  of  the  estimator  (23),  compar¬ 
ing  the  usage  of  the  quadratic  approximation  (25)  to  the  Poisson  estimator  (27).  Figures 
8  and  9  are  for  the  case  where  TV  =  7,  M  =  3,  l  =  2,  £o  =  0.001  and  r  =  1.  The  jth 
simulation  uses  K  =  10J .  Figure  8  compares  an  estimator  using  the  quadratic  approxima¬ 
tion  (25)  to  one  using  the  Poisson  estimator  (27),  with  H  =  1000.  In  this  example,  both 
estimators  seem  to  be  settling  down  after  K  =  1000  simulations.  Figure  9  is  for  the  same 
scenario,  except  only  estimators  based  upon  the  Poisson  estimate  (27)  are  included.  The 
three  cases  considered  are  for  H  =  K,  H  =  10  and  H  =  1000.  This  simulation,  as  well  as 
others  investigated,  showed  that  it  is  sufficient  to  take  around  H  =  1000  to  obtain  a  good 
Poisson  estimate. 

Figure  10  is  a  simulation  for  the  case  where  TV  =  5,  M  =  3,  l  =  2,  £o  =  1  and  r  =  5.3704. 
As  previously,  K  —  10-7  for  the  jth  simulation.  Three  estimators  are  compared  to  an 
approximation  based  upon  a  numerical  integration  scheme1.  The  numerical  integration 
scheme  gave  a  detection  probability  of  0.0095.  The  two  Poisson  estimators  use  H  =  K 
and  H  =  1000  respectively.  As  can  be  observed,  the  Poisson  based  estimates  coincide  with 
the  numerical  value  rather  quickly,  while  the  quadratic  based  estimator  improves  slowly. 

1This  was  provided  by  Mr  Daniel  Finch,  EWRD,  who  used  Matlab’ s  numerical  integration  function 
quad  to  estimate  the  detection  integral  (20). 
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Figure  11  contains  estimates  for  the  case  where  N  =  20,  M  =  15,  l  =  2,  £o  =  1  and 
t  =  1.4849.  The  three  estimators  used  are  the  same  as  in  Figure  10.  The  numerical 
scheme  gave  a  value  of  0.1184.  In  this  case,  the  Poisson  estimator  with  H  =  K  has  the 
best  performance. 

Other  simulations  considered  showed  that  the  estimator  (23),  when  coupled  with  the 
quadratic  approximation  (25),  had  very  poor  performance  when  the  threshold  r  and  the 
SNR  £0  were  fairly  large.  In  contrast  to  this,  it  was  found  that  using  the  Poisson  estimate 
(27)  improved  the  estimation  considerably,  with  reasonable  results  for  K  =  1000. 
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5  Conclusions  and  Future  Directions 


This  report  examined  the  Monte  Carlo  estimation  of  definite  integrals.  A  number  of 
estimators  of  the  Gamma  and  Beta  function  were  considered.  These  estimators  performed 
reasonably  well  in  practice.  A  number  of  estimators  of  the  probability  of  detection,  for  a 
binary  integration  scheme,  were  also  considered.  These  gave  reasonable  results  in  practice. 

In  future  work,  the  detection  probability  estimator  will  be  compared  to  other  estimators, 
to  gauge  its  performance. 
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Appendix  A:  Simulations 
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Estimators  for  Gamma  function  using  10000  samples 
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Figure  1:  Simulations  of  the  Gamma  estimators  (12)  and  (13),  for  a  selection  of  values 
of  X.  In  each  case,  both  estimators  are  compared  with  an  estimate  based  on  the  Matlab 
inbuilt  Gamma  function.  The  number  of  simulations  used  in  each  case  is  10,000.  In  the 
legend,  Gamma  refers  to  the  exact  result,  Exponential  refers  to  the  estimator  (12)  and 
Weibull  refers  to  (13). 
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Estimators  for  Gamma  function  using  10000  samples 
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Figure  2:  A  simulation  under  exactly  the  same  conditions  as  that  of  Figure  1 ,  showing 
slightly  worse  results. 
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Estimators  for  Gamma  function  using  100000  samples 
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Figure  3:  Simulations  of  Gamma  function  estimators  as  in  Figures  1  and 
number  of  simulation  samples  has  been  increased  to  100,000 . 
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Figure  4:  Simulations  of  the  Beta  integral  estimators  (If),  (15)  and  (16).  In  each  case, 
50  simulations  have  been  used  to  estimate  the  integral  (11)  with  parameters  (n,m).  In  the 
legend,  Beta  refers  to  the  exact  result,  while  the  other  three  refer  to  the  estimators  in  order 
of  appearance. 
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Figure  5:  Same  as  for  Figure  except  100  simulations  are  used  to  generate  the  estimates. 
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Estimators  for  Beta  function  using  1000  samples 
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Figure  6:  As  for  Figure  5 ,  except  1000  simulations  are  used. 
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Figure  1:  Simulation  of  (23)  in  the  case  of  N  =  5,  M  =  3,  1  =  2,  £o  =  1,  r  =  0A  and  for 
simulation  j,  K  =  10J.  The  quadratic  approximation  (25)  was  used  to  estimate  the  single 
pulse  probability  (19). 
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Figure  8:  A  comparison  of  estimator  (23)  using  the  quadratic  approximation  (25)  and  the 
Poisson  estimator  (27).  In  this  case,  N  —  7,  M  =  3;  l  =  2,  £o  —  0.001 ;  r  —  1.  For  each 
j,  K  =  10J  ;  and  77  =  1000. 
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Estimation  of  Detection  Probability 


Figure  9:  Comparison  of  three  estimators  using  different  Poissons  estimators.  The  same 
parameter  values  are  used  as  for  the  simulation  of  Figure  9 ,  except  the  value  of  H  is  as 
shown. 
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Figure  10:  A  simulation  in  the  case  where  N  =  5;  M  =  3;  l  =  2,  £o  =  1  and  T  —  5.3704. 
As  in  previous  simulations ,  K  —  10-7.  The  first  Poisson  estimate  uses  H  =  K ,  while 
the  second  uses  H  =  1000.  The  numerical  estimate  is  based  upon  numerical  integration 
applied  directly  to  the  detection  probability  integral  (20). 
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Figure  11:  As  for  Figure  10,  except  N  —  20,  M  =  15,  l  =  2,  £o  —  1  and  T  —  1.4849. 
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Appendix  B:  Matlab  Code 


Figure  12:  This  function  calculates  Matlab ’s  Gamma  function  (Option  0 ),  estimator  (12) 
(Option  1)  and  estimator  (13)  (Option  2).  N  Gamma  values  are  calculated  for  integers 
and  half-integers  from  1  to  N.  Parameters  a  and  b  correspond  to  the  a  and  (3  Gamma 
parameters.  Samples  is  a  vector  of  uniformly  sampled  points  from  the  unit  interval. 
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Figure  13:  Function  for  calculation  of  Matlab’s  Beta  function  (Option  0 ),  estimator  (If) 
(Option  1 ),  estimator  (15)  (Option  2)  and  estimator  (16)  (Option  3).  N  is  an  integer  indi¬ 
cating  the  size  of  the  matrix  of  Beta  values.  Samples  is  as  for  the  Gamma  implementation 
in  Figure  12.  AsMatrix ,  an  optional  argument ,  provides  a  better  formatted  output. 
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Figure  If:  Implementation  of  the  binary  integration  detection  probability  estimator  (23). 
Option  enables  the  single  pulse  probability  of  detection  to  be  estimated  via  (25)  (Option 
0)  or  (27)  (Option  1).  Parameter  Kl  is  the  number  of  Monte  Carlo  simulations  for  the 
main  estimator  (23),  while  K2  is  for  the  Monte  Carlo  estimator  (27).  N,  M,  l  and  tau 
are  the  corresponding  parameters  in  the  binary  integration  scheme ,  while  xO  =  £o- 
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